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We study the phase structure of QCD in the T — jl plane using a histogram method and the 
reweighting technique by performing phase quenched simulations of two-flavor QCD with RG- 
improved gauge action and 0(a) improved Wilson quark action. Taking the effects of the complex 
phase of the quark determinant using the cumulant expansion method, we calculate the probability 
distribution function of plaquette and phase-quenched determinant as a function of T and fl. We 
discuss the order of the QCD phase transition consulting the shape of the probability distribution 
function. 
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1. Introduction 

In order to reveal the phase structure of QCD, it is indispensable to study QCD by first principle 
lattice simulations. The lattice simulations, however, have the notorious sign problem at non-zero 
quark chemical potential ji. Here, we report our study on the phase structure of finite-density 
QCD using a histogram method with the reweighting technique [1,2] performing phase-quenched 
simulations and taking the effects of the complex phase of the quark determinant by the cumulant 
expansion. Results of test simulations have been reported in the previous Lattice conference [3]. 
We found that the method is promising. We have extended the study by increasing the number of 
simulation points as well as the statistics of each simulation point. By increasing ji/T, we found 
a preliminary signal suggesting the on- set of the first-order transition line separating the hadronic 
and QGP phases. 



2. Histogram method 

We define the probability distribution function w{0) for an operator O by the histogram of 
the measurement result O on the gauge configurations. We further define an effective potential by 
V(0) = — lnw(O). Near a first-order phase transition, w(0) will show a double-peak shape on 
finite lattices when O is sensitive to the phases. We may thus detect first-order transition from a 
negative curvature between two local minima in V(0). 

For simplicity, let us consider the case of the degenerate Nf flavors. Generalization to non- 
degenerate cases is straightforward. To explore the phase structure at finite T and /i, we choose 
the generalized plaquette P = — 5 g /6j8A^ S i te , where S g is the gauge action and 7V S i te = A^ 3 x N t is the 
space-time volume, and the absolute value of the quark determinant, = TVf In |detM (ju) / detM(O) |, 
as the operators for w. We have suppressed the arguments k and j8 in F and M because we fix k 
and csw in this study. 1 

Decomposing the quark determinant as (detM(/i)) iVf = e l6 ^ |detM(/i)| iVf , the partition func- 
tion at finite \i can be written as 

z(/?,/0 l 



J @Ue i6 ^ \detM(n)\ Nl e 6j3N ^ p = J dPdFw (P,F;P,n) 



Z(J8,0) Z(j8,0) 

' dPdFe-vW'-W, (2.1) 



where 



w (P,F;I3,h) = J mS(P-P)S(F-P(n)) IdetM^r^ 6 ^ (2-2) 

is the phase-quenched distribution function normalized by the partition function at jX = 0, and 

/ 8 (P - P) 8 (F - F {ii ) ) | detM( \i ) | Ni e 6 ^ p 



e ) (P,F;n) 



J 9U8{P - P)8{F - F(n )) |detM(jU ) \ Nt e 6 P N ^ p 
(e ie ^8(P-P)8(F-F(n)))) {M 



[S(P-P)S(F-F(ji 



(2.3) 



(i3,M) 



'When csw depends on ji, additional contributions to the reweighting have to be taken into account [4]. In this 
study, we instead fix c$w at all j5. This defines another renormalization scheme. Physical properties including the phase 
structure will not be affected. 
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is the expectation value of the complex phase of the quark determinant in phase-quenched simu- 
lation with fixed P and F. The double bracket is for the expectation value in the phase-quenched 
simulation. The effective potential is then given by V = — lnwo — ln(e l °). Note that (e l °) does not 
depend on j8 since the factor e 6 P NsiteP cancels out between numerator and denominator in Eq. (2.3). 

We carry out phase-quenched simulation at (j3o,jUo) and calculate wo at (j6,/i) using the 
reweighting technique as 



w (P,F;P,n)=R(P,F;P,Po,n,m>) w (P,F; j8 , Mo), 



/?(P,F,j8,ft,jU,jUo)=e' 



Similarly, the phase factor is given by 



8(P-P)8(F-F) 



detAf(ju) 
detM(/xo) 



N f 



[S(P-P)8(F-F))) 



(2.4) 



(2.5) 



(ft,Mo) 



e W )(P,F;n) 



JO 



detM(ju) 
detM(juo) 



8(P-P)8(F-F) 



detM(ju) 
detM(jUo) 



Aft 



8(P-P)8(F-F] 



(2.6) 



(ft,Mo) 



where the j6o -dependence in the right hand side cancels out. From these equations, we find that, 
under a shift of j8 at ji = /lo, the slope d V / changes by a constant factor, while <9 2 V / dP 2 remains 
the same. 



3. Cumulant expansion for the complex phase of the quark determinant 



Evaluation of the phase factor is suffered from the sign problem. To mitigate the problem, we 
apply the cumulant expansion method to evaluate the phase factor [1,5]: 



e ie )(P,F;n) = exp 



(3.1) 



The odd-order cumulants are the source of the sign problem. A key observation is that, because 
the system has the symmetry under \x <-> — fi, the odd-order cumulants should vanish when the 
statistics is sufficiently high. Keeping the even-power cumulants only, the phase factor is real and 
positive. Thus, the sign problem is resolved if the cumulant expansion converges. 

In a previous study, we found that, when we define 6 by the Taylor expansion of A^fJmlndetM( i u) 
in terms of /I, then the resulting distribution of is well described by the Gaussian function at small 
ji [1, 5]. Note that is defined in the range (— °o,°o), and the conventional 6 modulus 2n is re- 
produced by taking the principal value. With a Gaussian distribution, the cumulant expansion is 
dominated by the leading term (0 2 ) c , and thus the expansion converges. In this study with phase- 
quenched simulations, we adopt the definition 



fh 1 / 1 
Jo 



d(lndetM(/x)) 

Wm 



J A 



(3.2) 



i.e., on each configuration generated by the phase-quenched simulation at (j3o,jUo), we measure 
[<9lndetM/<9(jU/r)]£ for several values of jl, and calculate 0{fi) by integrating it numerically. For 
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Figure 1: The inverse susceptibilities of P (left panel) and F (right panel) as a function of j3 at various jtto/7 '. 

a smooth interpolation in the numerical integration procedure, we combine the information of the 
second derivative [<9 2 lndetM/<9(ju/r) 2 ] too. We confirm that the resulting distribution of 9 is 
close to Gaussian at the simulation points we study [3]. 

Using the results of [dlndetM/d(ji/T)}^ and [(9 2 lndetM/(9(^/r) 2 ] we also calculate 
F(n) and the reweighting factor |detM(/i)/detM(jUo)| in Eqs. (2.5) and (2.6) as follows: 



F(n) = Nf\n 
C(jU) = iVfln 



detM(ju) 



detM(0) 
detM(ju) 



detM(^o) 



= Nf <He 
Jo 

=Nf mc 

JfJo/T 



d(lndetM(ju))" 

Wm _ 

a(lndetM(ju)) 
dOi/T) 



d 



(3.3) 
(3.4) 



To the lowest order of the cumulant expansion, the curvatures of the effective potential V in P 
and F directions are now given by 



d 2 V 
dP 2 
d 2 V 
dF 2 



d 2 (— In wo) 
dP 2 



(P,F;j8o,Mo) 



d 2 lnR 

~dP r 



1«9 2 (0 2X 



(P,F;P,h,n,Ho) + 



(P,F;P,n) = d -^^ (P,F^^ 



2 dP 2 

\d 2 (e 2 ) 



dF 2 



2 dF 2 



(3.5) 



If the curvature of V is always positive, the transition between confined and deconfined phases is 
not first order. On the other hand, if negative curvature appears in the critical region, it indicates 
on-set of first-order transition. 



4. Numerical simulations and the results 

We use the RG-improved Iwasaki gauge action and the 0(a) -improved Wilson quark action. 
We generated gauge configurations on 8 3 x 4 lattice in two-flavor QCD with the phase-quenched 
measure at various values of (j3o, jUo) in the range j3o = 1.2-2.0 and (Jq/T = 2.0-4.0. K and csw are 
fixed to 0.141139 and 1.603830 respectively. We measure the first and the second derivatives of 
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Figure 2: Leftl: the statistical average of the second order cumulant of the complex phase as a function of 
p. Right: Schematic summary of the situation in the T — jl plane. Red dots represent the simulation points. 



IndetM every 10 trajectories using the random noise method of [5] with 50 noises. The statistics is 
2900 x 10 trajectories after thermalization. 

We find that the phase-quenched distributions for P and F around (P)^ jUo ) and (F)^ Q ^ have 
simple single bump forms. By approximating them by Gaussian functions, we estimate the second 
derivatives of In wo in Eq. (3.5) by inverse susceptibilities of P and F as 

d 2 (-lnwo) = 1 = 1 d 2 (-\nw ) = 1 1 

dP 2 ~ <(P-(P)) 2 ) xi dF 2 ~ ((F-(F)) 2 ) X i ( " } 

where ( • ) is the phase-quenched average at (/3o,jUo). 

The inverse susceptibilities are shown in Fig. 1 as function j8 at various /lo- We find that 
Xf increases rapidly as /lo increases. At small j6, ((F — (F)) 2 ) at /Jq/T = 3.6 is about 10 times 
smaller than that at jio/T = 2.4. A model calculation with an isotriplet chemical potential, which 
corresponds to the phase quenched two-flavor QCD, suggests that there is a second-order phase 
transition to the pion condensed phase at large /lo [6]. The increase of Xf suggests that the system 
is getting closer to the pion condensed phase as jio is increased. On the other hand, Xp is insensitive 
to /io- It is hard to expect a region with negative curvature in the P direction. We thus concentrate 
on the curvature in the F direction in the followings. 

Fig. 3 shows the results of d 2 (—lnwo)/dF 2 evaluated at various (/3o,/io) in the P-F plane 
for ji/T = 3.2, 3.6, and 4.0. The distribution function gives positive contribution and decreases 
quickly as the chemical potential increases, in the range we have studied. 

To evaluate InR and (0 2 ) c with Eqs. (2.5) and (2.6), we calculate the reweighting factor 
|detM(jii)/ detM(/io)| iVf = e c . We note that, in the range of parameters we have studied except 
just at ji = /lo where the reweighting factor is not needed, C shows a strong linear correlation 
with F such that C w aF + b is well satisfied in the ensemble at (j3o,jUo). When we approxi- 
mate C « aF + b, the reweighting factor e c can be factored out from the ensemble average with 
constrained F. Then, the resulting InR is a linear function of F, and InR does not contribute to 
d 2 V/dF 2 . 



5 



Phase structure of finite density QCD with a histogram method 



Y. Nakagawa 



2.0 r 



*=0.144139, n/T=3.2 



k=0.144139, fi/T=3.6 



K=0A44\39,jd/T=4.0 





100 200 300 400 500 600 700 ' 200 400 600 




0.00075 



1000 200 400 600 800 1000 1200 1400 



Figure 3: Contour plot of the curvature of the histogram, d 2 ( — lnwo)/dF 2 , at jl/T = 3.2,3.6,4.0. 
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Figure 4: Contour plot of the second order cumulant \ (0 2 ) c at ji/T = 3.2,3.6,4.0. 



With this approximation, we can also evaluate (0 2 ) c without the reweighting factor. In this 
study, we roughly estimate (0 2 ) c by the statistical average at the simulation point. In the study of 
curvatures, this corresponds to an approximation that 9 depends on P and F at most linearly in the 
dominant range around (P) and (F) at each simulation point. The results of (0 2 ) c are plotted in 
the left panel of Fig. 2 as a function of j8. We observe that (0 2 ) c shows a sharp bend downwards 
near the pseudo critical line when /io is large. When we translate the j8 -dependence to F and P 
dependences, this means that (0 2 ) c gives a negative contribution to the curvature of V. The second 
order cumulant in P — F plane is drawn in Fig. 4. We find that the fluctuation of the phase of the 
quark determinant increases in the small region, and the contours get closer to each other at 
larger fi/T. When the negative contribution becomes sufficiently strong, we may expect a negative 
curvature of V. This situation is depicted schematically in the right panel of Fig. 2. 

Fig. 5 shows the final results for the curvature of the effective potential in the F direction at 
jx/T = 3.2, 3.6, 4.0. We observe that the curvature remains positive in P — F plane at \ijT = 3.2 and 
3.6. At jx/T = 4.0, however, we see a region at small in which the curvature is consistent 

with zero. Although we could not find a region with negative curvature, the appearance of zero 
curvature region suggests that we have the end point of the first-order transition line separating the 
hadronic and QGP phases around this point. 
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Figure 5: Contour plot of the curvature of the effective potential in F direction at/i/r = 3.2,3.6,4.0. 



5. Summary 

We proposed a new approach to explore the phase structure of finite density QCD based on the 
histogram method and the reweighting technique using phase quenched simulations. We defined 
an effective potential in terms of the probability distribution function for a generalized plaquette 
(corresponding to the gauge action) and the absolute value of the quark determinant (corresponding 
to the quark action), and studied the curvature of the effective potential to detect the first order 
phase transition. We carried out a phase-quenched simulation of two-flavor QCD with improved 
Wilson quarks, and incorporated the effects of the complex phase of the quark determinant by 
the reweighting technique and the cumulant expansion method. We found that the complex phase 
factor contributes negatively to the curvature of the effective potential. If the curvature of the 
effective potential becomes negative around the minimum of the potential, the transition is first 
order. Although we could not find a region with negative curvature within the present range of 
simulation parameters, when we increase ji/T up to 4.0, we observe a region in which the curvature 
is consistent with zero. This suggests the end point of the first-order transition line separating the 
hadronic and QGP phases around this point. 
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